SCN5A Penetrance Prediction
Brett Kroncke and Yi Zuo
February 15, 2020
- Data preparation
- EM algorithm
- Step 1 Calculate weighted penetrance
- Step 2 Fit Beta-Binomial model using method of moments
- Step 3 Get posterior \(\alpha_{1i}\) and \(\beta_{1i}\)
- Step 4 Get posterior mean \(\tilde\mu_i\) and variance \(\tilde\sigma^2_i\)
- Step 5 Calculate penetrance density from \(\tilde\mu_i\) using function
funcdist - Step 6 Fit posterior mean \(\tilde\mu_i\) with clinical covariates
- Step 7 Use fitted penetrance and MSE from Step 6 to solve for \(\alpha\) and \(\beta\)
- Step 8 Calculate convergence criterion \(\delta\) to see if it’s small enough to end the loop
- Results
Data preparation
Data are cleaned and saved in RData format for the subsequent analyses.
# Step 1: load the original data
con = dbConnect(SQLite(),
dbname="VariantSCN5A-second-revision.db")
alltables = dbListTables(con)
my.data <- dbReadTable(con, 'VariantSCN5A')
my.data[my.data=='NA'] <- NA
data<-my.data
dbDisconnect(con)
dim(data) # 2417 by 29
names(data)
# clean `resnum`
table(data$resnum)
sum(is.na(data$resnum)) # 51 * + 1246 missing = 1297 abnormal
sum(data$resnum=="") # 1246 missing
# convert strings to numeric values
data$resnum<-suppressWarnings(as.integer(data$resnum))
sum(is.na(data$resnum)) # 1300 missing, three more
# replace missing resnum with the number after the first capital letters?
for (variant in data$var){
if (is.na(data[data$var==variant, "resnum"])){
suppressWarnings(data[data$var==variant, "resnum"]<-as.integer(strsplit(variant,"[A-Z]")[[1]][2]))
}
}
sum(is.na(data$resnum)) # 45 unsuccessful conversion
data[is.na(data$resnum),c("var","resnum")]
# manually get resnum from var among variants of interest
data[data$var=="E444fsX14","resnum"]<-444
data[data$var=="F1775fsX1785","resnum"]<-1775
data[data$var=="L1579fsX53","resnum"]<-1579
data[data$var=="p.1795insD","resnum"]<-1795
data[data$var=="p.G1031fsX27","resnum"]<-1031
data[data$var=="p.G1748del","resnum"]<-1748
data[data$var=="p.G1758del","resnum"]<-1758
data[data$var=="p.L1339del","resnum"]<-1339
data[data$var=="p.M741_T742delinsI ","resnum"]<-741
data[data$var=="Y1795insD","resnum"]<-1795
data[data$var=="p.A586_L587del","resnum"]<-586
data[data$var=="p.E1939_E1943del","resnum"]<-1939
data[data$var=="p.E1072del","resnum"]<-1072
data[data$var=="p.E1064del","resnum"]<-1064
data[data$var=="p.F1465_L1480dup","resnum"]<-1465
data[data$var=="p.M1875dup","resnum"]<-1875
data[data$var=="p.N507_L515dup","resnum"]<-507
data[data$var=="p.N291_S293dup","resnum"]<-291
data[data$var=="p.Pro2006del","resnum"]<-2006
data[data$var=="p.Gln1000del","resnum"]<-1000
data[data$var=="p.S1970_S1972del","resnum"]<-1970
data[data$var=="p.T290_G292del","resnum"]<-290
data[data$var=="p.I1758del","resnum"]<-1758
data[is.na(data$resnum),c("var","resnum")]
# after the conversion, how many missing in resnum?
sum(is.na(data$resnum)) # 23
# are resnum unique? No
length(unique(data$resnum))
# are nativeAA and mutAA different when resnum is the same?
duplicates <- data[duplicated(data$resnum)|duplicated(data$resnum,fromLast = T),c("var","resnum","nativeAA","mutAA")]
# it appears so
head(duplicates[order(duplicates$resnum),])
# nativeAA
table(data$nativeAA)
sum(is.na(data$nativeAA)) # 168 missing
# mutAA
table(data$mutAA)
sum(is.na(data$mutAA)) # 168 missing
# create a composite index for merging
data$res_nm <- paste0(as.character(data$resnum),data$nativeAA,data$mutAA)
# Step 2: load five supplementary data
# load data
bcl <- read.csv("covariates/scn5a_bcl_features.txt", header = TRUE, sep = "\t")
sasa <- read.csv("covariates/5a_5x0m_sasa.csv", header = TRUE)
pph2<-read.csv("covariates/pph2-20180329.csv", strip.white = TRUE)
prov<-read.csv("covariates/provean_revised.csv", strip.white = TRUE)
pam <- read.csv("covariates/sin_pam_energy_funcdist.csv", header = TRUE)
# Step 3: examine supplementary data
# pam
# pam has the same covariates as in the original data
sum(names(pam)%in%names(data)) == length(names(pam))
# seven times duplicated observations
nrow(pam) / nrow(unique(pam))
pam.u <- unique(pam)
head(pam.u)
data[data$var == "A2T","pamscore"]
# no need to merge with pam
# pph2
dim(pph2)
# all unique
dim(unique(pph2))
# contains four more covariates but no `var`
sum(names(pph2)%in%names(data)) == length(names(pph2))
names(pph2)%in%names(data)
names(pph2)
head(pph2)
# each residue number has 19 rows
length(unique(pph2$resnum))
pph2[pph2$resnum==1,]
table(data$resnum)
# create a composite index for merging
pph2$res_nm <- paste0(as.character(pph2$resnum),pph2$nativeAA,pph2$mutAA)
# sasa
dim(sasa)
# all unique
dim(unique(sasa))
head(sasa)
# unique resnum
length(unique(sasa$resnum)) == nrow(sasa)
# bcl
dim(bcl)
# eight duplicates
dim(unique(bcl))
head(bcl)
bcl[duplicated(bcl),]
bcl <- distinct(bcl)
dim(bcl)
# no duplicates
dim(unique(bcl))
head(bcl)
# prov
dim(prov)
dim(unique(prov))
# 15 duplicates
prov <- distinct(prov)
dim(unique(prov))[1] == nrow(prov)
head(prov)
# create a composite index for merging
prov$res_nm <- paste0(as.character(prov$resnum),prov$nativeAA,prov$mutAA)
# Step 4: merge data with 4 additional data sets
# pph
data.pph <- merge(data,pph2,by="res_nm",all.x=T)
dim(data.pph)
# drop the resnum, mutAA, nativeAA from pph2
data.pph<-data.pph[,-which(names(data.pph) %in% c("resnum.y","mutAA.y","nativeAA.y"))]
# rename resnum
colnames(data.pph)[colnames(data.pph)=="resnum.x"] <- "resnum"
colnames(data.pph)[colnames(data.pph)=="mutAA.x"] <- "mutAA"
colnames(data.pph)[colnames(data.pph)=="nativeAA.x"] <- "nativeAA"
names(data.pph)
# sasa
data.pph.sasa <- merge(data.pph,sasa,by="resnum",all.x=T)
dim(data.pph.sasa)
colnames(data.pph.sasa)
# bcl
data.pph.sasa.bcl <- merge(data.pph.sasa,bcl,by="var",all.x=T)
dim(data.pph.sasa.bcl)
names(data.pph.sasa.bcl)
# prov
data.pph.sasa.bcl.prov <- merge(data.pph.sasa.bcl,prov,by="res_nm",all.x = T)
dim(data.pph.sasa.bcl.prov)
names(data.pph.sasa.bcl.prov)
# drop the resnum, mutAA, nativeAA from prov
data.pph.sasa.bcl.prov<-data.pph.sasa.bcl.prov[,-which(names(data.pph.sasa.bcl.prov) %in% c("resnum.y","mutAA.y","nativeAA.y"))]
# rename resnum
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="resnum.x"] <- "resnum"
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="mutAA.x"] <- "mutAA"
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="nativeAA.x"] <- "nativeAA"
names(data.pph.sasa.bcl.prov)
dim(data.pph.sasa.bcl.prov)
head(data.pph.sasa.bcl.prov)
mydata <- data.pph.sasa.bcl.prov
# Step 5: exclusion criteria
# drop zero total carriers
mydata$gnomAD<-as.integer(mydata$gnomAD)
sum(is.na(mydata$gnomAD))
mydata$gnomAD[is.na(mydata$gnomAD)] <- 0
mydata$total_carriers<-mydata$unaff+mydata$lqt3+mydata$brs1+mydata$gnomAD
mydata.total<-mydata[mydata$total_carriers>0, ]
dim(mydata.total)
names(mydata.total)
# keep only missense, aadel and aains mututation types
mydata.total <- mydata.total[mydata.total$mut_type %in% c("missense","aadel", "aains"),]
dim(mydata.total)
# rename the final data set without clinical covariates
d <- mydata.total
d<-d[,!(names(d) %in% "aaneigh_rank_inc_loops")]
d<-d[!is.na(d$resnum),]
dim(d)
names(d)
save(d,file="BrS1_data.RData")EM algorithm
Step 1 Calculate weighted penetrance
The weighted penetrance is
\[\bar p_w=\frac{\sum_{i=1}^m w_i \hat p_i}{\sum_{i=1}^m w_i} \]
where the weight is
\[w_i=1-\frac{1}{0.01+n_i} \]
load("BrS1_data.RData")
# Step 1: Calculate overall penetrance and weighted penetrance
# observed penetrance for each variant
d$p_observed <- d$brs1/d$total_carriers
# overall penetrance
p_bar <- sum(d$brs1)/sum(d$total_carriers)
# weighted penetrance
wp <- function(p_observed,weight){
sum(weight*p_observed)/sum(weight)
}
# weighted penetrance
w <- 1-1/(0.01+d$total_carriers)
p_w <- wp(d$p_observed, w)Step 2 Fit Beta-Binomial model using method of moments
Variants have common prior
\[\hat p_i\sim Beta(\alpha_0,\beta_0) \]
The mean and variance of Beta distribution is
\[E(\hat p_i)=\bar p_w,\: Var(\hat p_i)=\text{MSE from weighted regression} \] We can solve for \(\alpha_0\) and \(\beta_0\) using method of moments.
# solve for alpha and beta in Beta distribution
solab <- function(mean, variance){
alpha <- (mean^2 * (1-mean) - variance * mean)/variance
beta <- alpha * (1 / mean - 1)
return(c(alpha,beta))
}
d$w <- w
# variance for weighted penetrance
vwp <- function(weight){
lm <- lm(p_observed~1,data=d,weights=weight)
mean((lm$residuals)^2*(lm$weights))
}
# weighted using weight 2
alpha_0_w <- solab(p_w, vwp(w))[1]
beta_0_w <- solab(p_w, vwp(w))[2]
# create empirical penetrance from alpha_0 and beta_0
d$p_empirical <- (d$brs1 + alpha_0_w)/(d$total_carriers+alpha_0_w+beta_0_w)Step 3 Get posterior \(\alpha_{1i}\) and \(\beta_{1i}\)
Variant-specific Beta-Binomial posterior \(\alpha_{1i}\) and \(\beta_{1i}\) based on Beta prior
\[\alpha_{1i}=\alpha_0+x_i,\beta_{1i}=\beta_0+n_i-x_i \]
Step 4 Get posterior mean \(\tilde\mu_i\) and variance \(\tilde\sigma^2_i\)
The mean \(\tilde\mu_i\) of posterior distribution is
\[\tilde \mu_i=\frac{\alpha_{1i}}{\alpha_{1i}+\beta_{1i}} \]
The variance \(\tilde\sigma^2_i\) of posterior distribution is
\[\tilde\sigma^2_i=\frac{\alpha_{1i}\beta_{1i}}{(\alpha_{1i}+\beta_{1i})^2(\alpha_{1i}+\beta_{1i}+1)} \]
which is raised at 20th percentile.
pmv <- function(alpha,beta){
mean <- alpha/(alpha+beta)
variance <- alpha*beta/(alpha+beta)^2/(alpha+beta+1)
return(cbind(mean,variance))
}
p_mean_w <- pmv(alpha_1_w,beta_1_w)[,1]
p_variance_w <- pmv(alpha_1_w,beta_1_w)[,2]
p_variance_w[p_variance_w<quantile(p_variance_w,0.2)] <- quantile(p_variance_w,0.2)
d$p_variance_w <- p_variance_wStep 5 Calculate penetrance density from \(\tilde\mu_i\) using function funcdist
Note that only 774 out of 1439 (53.8%) variants fall within structured regions of the protein.
index <- seq(1,length(d$resnum),1)
d$p_mean_w <- p_mean_w
w_results <- sapply(index,function(x) funcdist(d[x, "resnum"], d[x, "var"], d, distances, "p_mean_w", "sigmoid", 7))## [1] "p.Gln1000del"
## [1] "Q1000L"
## [1] "P1002S"
## [1] "C1004R"
## [1] "I1005T"
## [1] "I1005V"
## [1] "A1006S"
## [1] "T1007I"
## [1] "T1007N"
## [1] "P1008S"
## [1] "P1011L"
## [1] "P1011S"
## [1] "P1014S"
## [1] "E1015K"
## [1] "T1016M"
## [1] "K1018E"
## [1] "T101I"
## [1] "P1021S"
## [1] "R1023C"
## [1] "R1023H"
## [1] "R1023P"
## [1] "K1024R"
## [1] "E1025A"
## [1] "R1027P"
## [1] "R1027Q"
## [1] "R1027W"
## [1] "E1029K"
## [1] "E1032D"
## [1] "E1032K"
## [1] "Q1033R"
## [1] "P1034T"
## [1] "G1035V"
## [1] "G1040R"
## [1] "D1041G"
## [1] "D1041N"
## [1] "E1043K"
## [1] "V1045M"
## [1] "P1048S"
## [1] "R104G"
## [1] "R104Q"
## [1] "R104W"
## [1] "A1050T"
## [1] "V1051A"
## [1] "A1052D"
## [1] "E1053K"
## [1] "D1055G"
## [1] "T1056A"
## [1] "E1061D"
## [1] "D1062H"
## [1] "E1063G"
## [1] "p.E1064del"
## [1] "S1066G"
## [1] "L1067R"
## [1] "G1068A"
## [1] "G1068D"
## [1] "T1069M"
## [1] "S106T"
## [1] "E1071K"
## [1] "p.E1072del"
## [1] "S1074G"
## [1] "S1074R"
## [1] "S1079F"
## [1] "S1079T"
## [1] "S1079Y"
## [1] "V1082A"
## [1] "S1083C"
## [1] "G1084D"
## [1] "G1084R"
## [1] "G1084S"
## [1] "A1088T"
## [1] "A1088V"
## [1] "P1090L"
## [1] "P1090Q"
## [1] "D1091A"
## [1] "D1091Y"
## [1] "W1095C"
## [1] "S1096C"
## [1] "S1096G"
## [1] "Q1097H"
## [1] "V1098L"
## [1] "V1098M"
## [1] "N109K"
## [1] "T10S"
## [1] "A1100T"
## [1] "A1100V"
## [1] "A1102T"
## [1] "S1103F"
## [1] "S1103Y"
## [1] "E1105V"
## [1] "A1106T"
## [1] "E1107K"
## [1] "S1109G"
## [1] "A110T"
## [1] "A1113T"
## [1] "A1113V"
## [1] "D1114E"
## [1] "D1114N"
## [1] "W1115R"
## [1] "R1116Q"
## [1] "R1116W"
## [1] "A1121V"
## [1] "A1125G"
## [1] "A1125T"
## [1] "A1125V"
## [1] "G1129S"
## [1] "Y112C"
## [1] "E1130K"
## [1] "T1131I"
## [1] "T1131S"
## [1] "P1132S"
## [1] "D1134E"
## [1] "D1134N"
## [1] "S1135I"
## [1] "C1136Y"
## [1] "V113A"
## [1] "V113I"
## [1] "S1140T"
## [1] "T1147I"
## [1] "T1147N"
## [1] "T1147S"
## [1] "A1148S"
## [1] "A1148T"
## [1] "Q1153H"
## [1] "I1154N"
## [1] "P1155S"
## [1] "D1156G"
## [1] "G1158S"
## [1] "S115G"
## [1] "D1163E"
## [1] "D1163G"
## [1] "D1163Y"
## [1] "P1164T"
## [1] "E1165D"
## [1] "E1165Q"
## [1] "D1166N"
## [1] "C1167Y"
## [1] "F1168L"
## [1] "T1169I"
## [1] "V1173D"
## [1] "R1174G"
## [1] "R1174W"
## [1] "R1175H"
## [1] "P1177L"
## [1] "C1178Y"
## [1] "A1180V"
## [1] "V1181A"
## [1] "V1181L"
## [1] "V1181M"
## [1] "T1183I"
## [1] "A1186T"
## [1] "P1187Q"
## [1] "K1189T"
## [1] "V1190F"
## [1] "R1193Q"
## [1] "R1193W"
## [1] "L1194M"
## [1] "R1195H"
## [1] "R1195S"
## [1] "Y1199S"
## [1] "P119L"
## [1] "P119S"
## [1] "S11R"
## [1] "H1200R"
## [1] "H1200Y"
## [1] "I1201M"
## [1] "V1202M"
## [1] "E1208K"
## [1] "T1209R"
## [1] "F1210S"
## [1] "p.I1212del"
## [1] "M1214T"
## [1] "I1215V"
## [1] "L1216V"
## [1] "S1218I"
## [1] "S1218T"
## [1] "S1219N"
## [1] "R121Q"
## [1] "R121W"
## [1] "G1220E"
## [1] "A1221V"
## [1] "L1222R"
## [1] "E1225K"
## [1] "G1225K"
## [1] "Y1228C"
## [1] "Y1228F"
## [1] "Y1228H"
## [1] "E1230K"
## [1] "E1231K"
## [1] "R1232Q"
## [1] "R1232W"
## [1] "K1233E"
## [1] "K1236N"
## [1] "K1236R"
## [1] "V1237F"
## [1] "L1239P"
## [1] "A123G"
## [1] "A123V"
## [1] "E1240Q"
## [1] "D1243N"
## [1] "K1244E"
## [1] "M1245I"
## [1] "T1247I"
## [1] "V1249D"
## [1] "A124D"
## [1] "V1251M"
## [1] "E1253G"
## [1] "L1255M"
## [1] "V125L"
## [1] "A1260D"
## [1] "G1262S"
## [1] "K1264N"
## [1] "K1264R"
## [1] "T1268N"
## [1] "T1268S"
## [1] "N1269S"
## [1] "K126E"
## [1] "A1270S"
## [1] "W1271C"
## [1] "W1273C"
## [1] "D1275N"
## [1] "I1278N"
## [1] "V1279I"
## [1] "V1281F"
## [1] "S1282A"
## [1] "L1283M"
## [1] "A1288G"
## [1] "F1293S"
## [1] "A1294G"
## [1] "E1295K"
## [1] "M1296T"
## [1] "P1298L"
## [1] "R1303Q"
## [1] "R1303W"
## [1] "T1304M"
## [1] "R1306H"
## [1] "R1306S"
## [1] "L1308F"
## [1] "R1309C"
## [1] "R1309H"
## [1] "L1311P"
## [1] "R1316L"
## [1] "R1316Q"
## [1] "F1317C"
## [1] "G1319V"
## [1] "M1320I"
## [1] "R1321K"
## [1] "V1323G"
## [1] "N1325S"
## [1] "A1326S"
## [1] "V1328M"
## [1] "G1329S"
## [1] "A1330D"
## [1] "A1330P"
## [1] "A1330T"
## [1] "I1331V"
## [1] "P1332L"
## [1] "P1332Q"
## [1] "I1334V"
## [1] "M1335R"
## [1] "L1338V"
## [1] "L1339F"
## [1] "p.L1339del"
## [1] "V1340I"
## [1] "F1344L"
## [1] "F1344S"
## [1] "W1345C"
## [1] "L1346I"
## [1] "L1346P"
## [1] "F1348L"
## [1] "N134S"
## [1] "I1350L"
## [1] "I1350T"
## [1] "M1351R"
## [1] "M1351V"
## [1] "V1353M"
## [1] "A1357V"
## [1] "G1358R"
## [1] "G1358W"
## [1] "K1359M"
## [1] "K1359N"
## [1] "M135V"
## [1] "F1360C"
## [1] "R1362S"
## [1] "C1363Y"
## [1] "I1364V"
## [1] "N1365S"
## [1] "Q1366H"
## [1] "Q1366R"
## [1] "G1369R"
## [1] "L136P"
## [1] "D1370A"
## [1] "D1370G"
## [1] "Y1375H"
## [1] "I1377M"
## [1] "V1378M"
## [1] "I137V"
## [1] "N1380K"
## [1] "p.N1380del"
## [1] "S1382I"
## [1] "C1384Y"
## [1] "L1387F"
## [1] "M138I"
## [1] "G1391R"
## [1] "V1398M"
## [1] "p.I137_C139dup"
## [1] "V1400I"
## [1] "V1405L"
## [1] "V1405M"
## [1] "G1406E"
## [1] "G1406R"
## [1] "G1408R"
## [1] "Y1409C"
## [1] "L1412F"
## [1] "Q1414H"
## [1] "A1416E"
## [1] "A1416G"
## [1] "K1419E"
## [1] "I141N"
## [1] "I141V"
## [1] "G1420D"
## [1] "G1420P"
## [1] "G1420R"
## [1] "G1420V"
## [1] "M1422R"
## [1] "D1423H"
## [1] "I1424V"
## [1] "A1427E"
## [1] "A1427S"
## [1] "A1428S"
## [1] "A1428V"
## [1] "D1430N"
## [1] "S1431C"
## [1] "R1432G"
## [1] "R1432S"
## [1] "G1433R"
## [1] "G1433V"
## [1] "G1433W"
## [1] "P1438L"
## [1] "Q1439H"
## [1] "Q1439R"
## [1] "E1441Q"
## [1] "Y1442C"
## [1] "Y1442N"
## [1] "N1443S"
## [1] "L1444I"
## [1] "Y1445H"
## [1] "I1448L"
## [1] "I1448T"
## [1] "Y1449C"
## [1] "Y1449S"
## [1] "V1451D"
## [1] "V1451L"
## [1] "S1458Y"
## [1] "F1460L"
## [1] "T1461S"
## [1] "N1463Y"
## [1] "L1464P"
## [1] "p.F1465_L1480dup"
## [1] "V1468A"
## [1] "V1468F"
## [1] "I1469V"
## [1] "V146A"
## [1] "V146M"
## [1] "N1472S"
## [1] "p.N1472del"
## [1] "F1473C"
## [1] "F1473S"
## [1] "Q1475L"
## [1] "Q1476R"
## [1] "K1477N"
## [1] "K1478E"
## [1] "G1481E"
## [1] "G1481R"
## [1] "G1481V"
## [1] "Q1483H"
## [1] "F1486L"
## [1] "p.F1486del"
## [1] "M1487K"
## [1] "M1487L"
## [1] "T1488R"
## [1] "E1489D"
## [1] "Q1491H"
## [1] "K1493R"
## [1] "p.K1493del"
## [1] "Y1495S"
## [1] "M1498R"
## [1] "M1498T"
## [1] "M1498V"
## [1] "R14C"
## [1] "R14H"
## [1] "R14S"
## [1] "p.K1500del"
## [1] "L1501V"
## [1] "p.L1501_K1505del"
## [1] "G1502A"
## [1] "G1502S"
## [1] "S1503Y"
## [1] "K1504E"
## [1] "K1505N"
## [1] "p.K1505_Q1507del"
## [1] "P1506S"
## [1] "P1506T"
## [1] "p.Q1507_P1509del"
## [1] "P1509T"
## [1] "R1512L"
## [1] "R1512Q"
## [1] "R1512W"
## [1] "L1514M"
## [1] "N1515S"
## [1] "I1521K"
## [1] "I1521T"
## [1] "D1523N"
## [1] "I1524T"
## [1] "V1525A"
## [1] "V1525M"
## [1] "T1526P"
## [1] "K1527R"
## [1] "D152N"
## [1] "V1532F"
## [1] "V1532I"
## [1] "T1533I"
## [1] "C1539F"
## [1] "C1539Y"
## [1] "V1543A"
## [1] "V1543L"
## [1] "T1544P"
## [1] "M1546T"
## [1] "V1547L"
## [1] "E1548K"
## [1] "G1548K"
## [1] "P154L"
## [1] "D1551N"
## [1] "D1551Y"
## [1] "Q1552L"
## [1] "Q1552R"
## [1] "S1553R"
## [1] "E1555K"
## [1] "I1557V"
## [1] "I1559V"
## [1] "L1560F"
## [1] "L1565M"
## [1] "F1567L"
## [1] "A1569P"
## [1] "W156R"
## [1] "I1570V"
## [1] "p.1570_F1571insI"
## [1] "p.I1570dup"
## [1] "F1571C"
## [1] "E1574K"
## [1] "C1575S"
## [1] "T157I"
## [1] "A1581S"
## [1] "L1582P"
## [1] "R1583C"
## [1] "R1583H"
## [1] "Y1585C"
## [1] "F1587V"
## [1] "T1588I"
## [1] "K158T"
## [1] "I1593M"
## [1] "F1594S"
## [1] "F1596C"
## [1] "F1596I"
## [1] "V1597M"
## [1] "V1598A"
## [1] "Y159C"
## [1] "R15G"
## [1] "R15M"
## [1] "R15T"
## [1] "L1601H"
## [1] "I1603F"
## [1] "V1604M"
## [1] "G1605C"
## [1] "G1605D"
## [1] "T1606I"
## [1] "S1609L"
## [1] "S1609W"
## [1] "D1610G"
## [1] "I1611V"
## [1] "Q1613H"
## [1] "Q1613L"
## [1] "p.F1617del"
## [1] "P1619L"
## [1] "P1619Q"
## [1] "E161K"
## [1] "E161Q"
## [1] "T1620K"
## [1] "T1620M"
## [1] "R1623L"
## [1] "R1623Q"
## [1] "V1624I"
## [1] "R1626C"
## [1] "R1626H"
## [1] "R1626L"
## [1] "R1626P"
## [1] "R1629G"
## [1] "R1629Q"
## [1] "Y162C"
## [1] "Y162H"
## [1] "I1630R"
## [1] "I1630V"
## [1] "G1631D"
## [1] "R1632C"
## [1] "R1632H"
## [1] "R1632L"
## [1] "L1634P"
## [1] "R1638Q"
## [1] "G1639A"
## [1] "G1642E"
## [1] "I1643L"
## [1] "R1644C"
## [1] "R1644H"
## [1] "R1644L"
## [1] "T1645M"
## [1] "A1649V"
## [1] "F164L"
## [1] "L1650F"
## [1] "M1652R"
## [1] "M1652T"
## [1] "I1660S"
## [1] "I1660V"
## [1] "G1661E"
## [1] "G1661R"
## [1] "V1667I"
## [1] "M1668T"
## [1] "A166T"
## [1] "S1672Y"
## [1] "F1674V"
## [1] "M1676I"
## [1] "M1676T"
## [1] "N1678S"
## [1] "A1680P"
## [1] "A1680T"
## [1] "Y1681F"
## [1] "W1684R"
## [1] "D1689N"
## [1] "D1690N"
## [1] "A1698T"
## [1] "F16L"
## [1] "M1701I"
## [1] "L1704H"
## [1] "Q1706H"
## [1] "T1708I"
## [1] "p.T1709del"
## [1] "T1709M"
## [1] "T1709R"
## [1] "F170I"
## [1] "S1710L"
## [1] "G1712C"
## [1] "G1712S"
## [1] "D1714G"
## [1] "L1717P"
## [1] "S1718R"
## [1] "N1722D"
## [1] "T1723N"
## [1] "P1725L"
## [1] "C1728R"
## [1] "C1728W"
## [1] "C1728Y"
## [1] "D1729N"
## [1] "P1730A"
## [1] "P1730H"
## [1] "P1730L"
## [1] "G1737D"
## [1] "S1738F"
## [1] "S1738T"
## [1] "R1739Q"
## [1] "R1739W"
## [1] "G1740R"
## [1] "D1741E"
## [1] "D1741N"
## [1] "D1741Y"
## [1] "G1743E"
## [1] "G1743R"
## [1] "S1744I"
## [1] "A1746T"
## [1] "A1746V"
## [1] "V1747M"
## [1] "G1748D"
## [1] "p.G1748del"
## [1] "I1749N"
## [1] "V174I"
## [1] "L1750F"
## [1] "T1753A"
## [1] "I1756V"
## [1] "I1758V"
## [1] "p.I1758del"
## [1] "S1759C"
## [1] "K175N"
## [1] "L1761F"
## [1] "L1761H"
## [1] "I1762M"
## [1] "p.I1762del"
## [1] "V1763L"
## [1] "V1763M"
## [1] "V1764F"
## [1] "M1766L"
## [1] "M1766T"
## [1] "M1766V"
## [1] "Y1767C"
## [1] "I1768V"
## [1] "I1770V"
## [1] "I1771T"
## [1] "L1772V"
## [1] "N1774D"
## [1] "F1775V"
## [1] "V1777L"
## [1] "V1777M"
## [1] "T1779M"
## [1] "L177P"
## [1] "E1780G"
## [1] "E1781D"
## [1] "E1781G"
## [1] "E1784K"
## [1] "L1786Q"
## [1] "L1786R"
## [1] "S1787N"
## [1] "A178G"
## [1] "D1790G"
## [1] "D1790N"
## [1] "p.D1790del"
## [1] "D1792N"
## [1] "D1792V"
## [1] "D1792Y"
## [1] "M1793K"
## [1] "p.Y1795_E1796insD"
## [1] "Y1795C"
## [1] "Y1795H"
## [1] "Y1795N"
## [1] "I1797V"
## [1] "R179Q"
## [1] "I1809M"
## [1] "G180V"
## [1] "Y1811N"
## [1] "S1812L"
## [1] "D1816E"
## [1] "D1816N"
## [1] "D1819N"
## [1] "A1820T"
## [1] "A1820V"
## [1] "E1823K"
## [1] "P1824A"
## [1] "L1825P"
## [1] "R1826C"
## [1] "R1826H"
## [1] "A1828S"
## [1] "A1828T"
## [1] "C182R"
## [1] "C182Y"
## [1] "Q1832E"
## [1] "I1833M"
## [1] "S1834R"
## [1] "L1835F"
## [1] "I1836T"
## [1] "D1839G"
## [1] "M1842L"
## [1] "M1842T"
## [1] "M1842V"
## [1] "G1845R"
## [1] "R1847C"
## [1] "R1847H"
## [1] "H1849R"
## [1] "H184R"
## [1] "C1850S"
## [1] "M1851I"
## [1] "M1851V"
## [1] "D1852V"
## [1] "I1853V"
## [1] "A185T"
## [1] "A185V"
## [1] "V1861F"
## [1] "V1861I"
## [1] "A1870T"
## [1] "K1872N"
## [1] "I1873V"
## [1] "M1875K"
## [1] "M1875T"
## [1] "p.M1875dup"
## [1] "E1877K"
## [1] "T187I"
## [1] "T187S"
## [1] "M1880V"
## [1] "P1884L"
## [1] "Y1889C"
## [1] "E1890K"
## [1] "T1895I"
## [1] "L1896P"
## [1] "R1897Q"
## [1] "R1897W"
## [1] "R1897Y"
## [1] "R1898C"
## [1] "R1898H"
## [1] "R18Q"
## [1] "R18W"
## [1] "E1901K"
## [1] "E1901Q"
## [1] "E1902A"
## [1] "V1903M"
## [1] "S1904L"
## [1] "M1906T"
## [1] "M1906V"
## [1] "I1908V"
## [1] "Q1909R"
## [1] "R190G"
## [1] "R190L"
## [1] "R190Q"
## [1] "R190W"
## [1] "R1910K"
## [1] "R1913C"
## [1] "R1913H"
## [1] "R1913S"
## [1] "R1914G"
## [1] "H1915N"
## [1] "H1915P"
## [1] "H1915Q"
## [1] "H1915Y"
## [1] "R1919C"
## [1] "R1919H"
## [1] "S1920C"
## [1] "K1922N"
## [1] "K1922R"
## [1] "H1923D"
## [1] "H1923Y"
## [1] "A1924T"
## [1] "S1925F"
## [1] "L1927P"
## [1] "F1928V"
## [1] "R1929C"
## [1] "R1929H"
## [1] "Q1930H"
## [1] "A1932V"
## [1] "G1933A"
## [1] "G1933D"
## [1] "G1933V"
## [1] "G1935S"
## [1] "S1937A"
## [1] "E1938K"
## [1] "p.E1939_E1943del"
## [1] "W193R"
## [1] "P1942H"
## [1] "P1942S"
## [1] "R1944Q"
## [1] "I1948V"
## [1] "A1949S"
## [1] "A1949T"
## [1] "Y1950C"
## [1] "V1951L"
## [1] "V1951M"
## [1] "E1954K"
## [1] "N1955Y"
## [1] "S1957P"
## [1] "R1958P"
## [1] "R1958Q"
## [1] "P1962L"
## [1] "P1962S"
## [1] "P1963L"
## [1] "S1964F"
## [1] "S1965G"
## [1] "S1965N"
## [1] "I1968M"
## [1] "I1968N"
## [1] "I1968S"
## [1] "I1968T"
## [1] "I1968V"
## [1] "p.S1970_S1972del"
## [1] "F1973L"
## [1] "P1975T"
## [1] "S1976C"
## [1] "Y1977N"
## [1] "V1980F"
## [1] "R1982T"
## [1] "A1983G"
## [1] "A1983V"
## [1] "T1984I"
## [1] "S1985R"
## [1] "D1986G"
## [1] "D1986N"
## [1] "N1987K"
## [1] "L1988R"
## [1] "V1990L"
## [1] "R1991Q"
## [1] "R1991W"
## [1] "G1992A"
## [1] "S1996N"
## [1] "S1996R"
## [1] "H1997R"
## [1] "S199T"
## [1] "M1I"
## [1] "D2000Y"
## [1] "A2002T"
## [1] "D2003N"
## [1] "F2004I"
## [1] "F2004L"
## [1] "F2004V"
## [1] "p.F2004dup"
## [1] "P2005A"
## [1] "P2005L"
## [1] "P2005S"
## [1] "p.Pro2006del"
## [1] "P2006A"
## [1] "P2006T"
## [1] "P2008L"
## [1] "D2009E"
## [1] "R2010G"
## [1] "R2012C"
## [1] "R2012H"
## [1] "V2016L"
## [1] "V2016M"
## [1] "I202T"
## [1] "A204T"
## [1] "A204V"
## [1] "E208K"
## [1] "N209S"
## [1] "N209T"
## [1] "S20F"
## [1] "I210T"
## [1] "L212P"
## [1] "L212Q"
## [1] "S216L"
## [1] "R219C"
## [1] "R219H"
## [1] "L21V"
## [1] "T220I"
## [1] "R222L"
## [1] "R222Q"
## [1] "V223L"
## [1] "L224F"
## [1] "R225Q"
## [1] "R225W"
## [1] "A226G"
## [1] "A226V"
## [1] "L227P"
## [1] "K228R"
## [1] "A22V"
## [1] "I230M"
## [1] "I230T"
## [1] "I230V"
## [1] "V232F"
## [1] "V232I"
## [1] "P234S"
## [1] "G235R"
## [1] "I239V"
## [1] "I239V "
## [1] "A23S"
## [1] "V240M"
## [1] "A242D"
## [1] "Q245K"
## [1] "V247L"
## [1] "V258A"
## [1] "E25K"
## [1] "S262G"
## [1] "V263I"
## [1] "A265V"
## [1] "L266H"
## [1] "G268S"
## [1] "Q270K"
## [1] "L271V"
## [1] "G274C"
## [1] "N275K"
## [1] "L276P"
## [1] "L276Q"
## [1] "H278D"
## [1] "H278R"
## [1] "R27C"
## [1] "R27H"
## [1] "R27L"
## [1] "C280Y"
## [1] "V281M"
## [1] "R282C"
## [1] "R282H"
## [1] "T285K"
## [1] "A286S"
## [1] "A286V"
## [1] "N288S"
## [1] "G289S"
## [1] "M28I"
## [1] "M28L"
## [1] "M28T"
## [1] "p.T290_G292del"
## [1] "p.N291_S293dup"
## [1] "N291H"
## [1] "N291K"
## [1] "N291S"
## [1] "G292S"
## [1] "V294M"
## [1] "E295K"
## [1] "G298D"
## [1] "G298S"
## [1] "L299F"
## [1] "L299M"
## [1] "A29E"
## [1] "A29V"
## [1] "A2T"
## [1] "V300I"
## [1] "D305N"
## [1] "L306F"
## [1] "L306V"
## [1] "S309C"
## [1] "S309N"
## [1] "E30G"
## [1] "D310N"
## [1] "E312K"
## [1] "L315P"
## [1] "K317E"
## [1] "K317M"
## [1] "K317N"
## [1] "G319C"
## [1] "G319R"
## [1] "G319S"
## [1] "T320N"
## [1] "S321Y"
## [1] "L325R"
## [1] "S330F"
## [1] "A332T"
## [1] "C335R"
## [1] "C335S"
## [1] "P336L"
## [1] "R340Q"
## [1] "R340W"
## [1] "C341Y"
## [1] "A344S"
## [1] "E346D"
## [1] "E346G"
## [1] "E346K"
## [1] "P348A"
## [1] "D349N"
## [1] "R34C"
## [1] "R34H"
## [1] "H350Q"
## [1] "G351C"
## [1] "G351D"
## [1] "G351S"
## [1] "G351V"
## [1] "Y352C"
## [1] "T353I"
## [1] "F355C"
## [1] "F355I"
## [1] "D356N"
## [1] "A359T"
## [1] "G35S"
## [1] "R367C"
## [1] "R367H"
## [1] "R367L"
## [1] "M369K"
## [1] "T370M"
## [1] "Q371E"
## [1] "W374G"
## [1] "R376C"
## [1] "R376H"
## [1] "T37A"
## [1] "S384T"
## [1] "A385T"
## [1] "G386E"
## [1] "G386R"
## [1] "I388S"
## [1] "Y389H"
## [1] "V396A"
## [1] "V396L"
## [1] "I397F"
## [1] "I397T"
## [1] "I397V"
## [1] "L39F"
## [1] "N3K"
## [1] "N3S"
## [1] "G400A"
## [1] "G400E"
## [1] "G400R"
## [1] "S401P"
## [1] "F402L"
## [1] "L404Q"
## [1] "L404V"
## [1] "N406K"
## [1] "N406S"
## [1] "L409P"
## [1] "L409V"
## [1] "A410V"
## [1] "V411M"
## [1] "V412D"
## [1] "A413E"
## [1] "A413T"
## [1] "M414V"
## [1] "A415T"
## [1] "Y416C"
## [1] "E418K"
## [1] "I424M"
## [1] "A425P"
## [1] "A425T"
## [1] "E428K"
## [1] "E429K"
## [1] "p.E429del"
## [1] "K430E"
## [1] "R433C"
## [1] "R433H"
## [1] "R433S"
## [1] "A437V"
## [1] "M438L"
## [1] "M438T"
## [1] "E439K"
## [1] "E439V"
## [1] "R43Q"
## [1] "L441F"
## [1] "H445D"
## [1] "H445Q"
## [1] "H445Y"
## [1] "E446K"
## [1] "A447G"
## [1] "A447S"
## [1] "T449A"
## [1] "Y449C"
## [1] "G452D"
## [1] "V453M"
## [1] "V456M"
## [1] "R458C"
## [1] "R458H"
## [1] "S459G"
## [1] "G45A"
## [1] "L461V"
## [1] "E462A"
## [1] "E462K"
## [1] "M463R"
## [1] "M463T"
## [1] "L466F"
## [1] "P468L"
## [1] "V469I"
## [1] "N470K"
## [1] "R474G"
## [1] "R474K"
## [1] "R475K"
## [1] "R475S"
## [1] "K480N"
## [1] "R481Q"
## [1] "R481W"
## [1] "M482I"
## [1] "T486A"
## [1] "T486S"
## [1] "E48K"
## [1] "G490A"
## [1] "G490E"
## [1] "E491G"
## [1] "R493K"
## [1] "K496M"
## [1] "K496N"
## [1] "S497C"
## [1] "E49K"
## [1] "F4V"
## [1] "E500K"
## [1] "D501G"
## [1] "P503S"
## [1] "R504T"
## [1] "A505E"
## [1] "M506K"
## [1] "p.N507_L515dup"
## [1] "T512I"
## [1] "R513C"
## [1] "R513H"
## [1] "R513P"
## [1] "G514C"
## [1] "R517S"
## [1] "S519F"
## [1] "A51V"
## [1] "M520R"
## [1] "M520V"
## [1] "K521E"
## [1] "P522S"
## [1] "R523C"
## [1] "R523H"
## [1] "R523S"
## [1] "S524Y"
## [1] "S525G"
## [1] "R526C"
## [1] "R526H"
## [1] "G527R"
## [1] "S528R"
## [1] "P52H"
## [1] "P52S"
## [1] "F530S"
## [1] "F530V"
## [1] "T531A"
## [1] "T531I"
## [1] "F532C"
## [1] "F532L"
## [1] "R533C"
## [1] "R533H"
## [1] "R533S"
## [1] "R535G"
## [1] "R535Q"
## [1] "D536H"
## [1] "G538D"
## [1] "R53Q"
## [1] "F543L"
## [1] "A551T"
## [1] "A551V"
## [1] "G552R"
## [1] "G552W"
## [1] "E553K"
## [1] "S554I"
## [1] "S554N"
## [1] "E555K"
## [1] "H557L"
## [1] "H557Q"
## [1] "H557Y"
## [1] "H558R"
## [1] "T559I"
## [1] "T559R"
## [1] "V563G"
## [1] "L567Q"
## [1] "R568C"
## [1] "R568H"
## [1] "R569Q"
## [1] "R569W"
## [1] "T570N"
## [1] "S571I"
## [1] "A572D"
## [1] "A572F"
## [1] "A572S"
## [1] "A572V"
## [1] "Q573E"
## [1] "Q573R"
## [1] "G574E"
## [1] "S577N"
## [1] "P578R"
## [1] "P578T"
## [1] "G579R"
## [1] "S581L"
## [1] "P583L"
## [1] "G584R"
## [1] "A586G"
## [1] "A586T"
## [1] "p.A586_L587del"
## [1] "p.586_587delAL"
## [1] "H588N"
## [1] "H588R"
## [1] "K590Q"
## [1] "N592K"
## [1] "N592S"
## [1] "D596G"
## [1] "C597G"
## [1] "C597Y"
## [1] "G599R"
## [1] "V601A"
## [1] "L604V"
## [1] "G605E"
## [1] "A606T"
## [1] "G607D"
## [1] "G607R"
## [1] "G607V"
## [1] "D608H"
## [1] "D608N"
## [1] "A60P"
## [1] "E610K"
## [1] "P614S"
## [1] "G615E"
## [1] "S616N"
## [1] "L618F"
## [1] "L619F"
## [1] "R620C"
## [1] "R620H"
## [1] "M623T"
## [1] "L624Q"
## [1] "E625D"
## [1] "E625Q"
## [1] "P627L"
## [1] "P628R"
## [1] "D629Y"
## [1] "K62T"
## [1] "T630M"
## [1] "T632M"
## [1] "S634L"
## [1] "S634W"
## [1] "E636K"
## [1] "P637L"
## [1] "G638D"
## [1] "G639A"
## [1] "G639R"
## [1] "K63N"
## [1] "P640A"
## [1] "P640L"
## [1] "P640S"
## [1] "A647D"
## [1] "A647S"
## [1] "A647V"
## [1] "P648L"
## [1] "C649R"
## [1] "C649Y"
## [1] "p.L64del"
## [1] "D651H"
## [1] "G652D"
## [1] "G652S"
## [1] "E654D"
## [1] "E654K"
## [1] "E654Q"
## [1] "E655K"
## [1] "P656L"
## [1] "A658V"
## [1] "R659Q"
## [1] "R659W"
## [1] "R661Q"
## [1] "R661W"
## [1] "A662S"
## [1] "S664G"
## [1] "A665S"
## [1] "A665T"
## [1] "V668I"
## [1] "D66N"
## [1] "S671C"
## [1] "A672S"
## [1] "A672T"
## [1] "L673P"
## [1] "L673V"
## [1] "E677V"
## [1] "R680C"
## [1] "R680H"
## [1] "H681P"
## [1] "C683G"
## [1] "C683R"
## [1] "C683S"
## [1] "C686F"
## [1] "R689C"
## [1] "R689H"
## [1] "Y68C"
## [1] "A691S"
## [1] "A691T"
## [1] "Q692K"
## [1] "R693C"
## [1] "R693H"
## [1] "Y694C"
## [1] "I696N"
## [1] "G69D"
## [1] "L6S"
## [1] "P701L"
## [1] "M704T"
## [1] "S705F"
## [1] "S705P"
## [1] "G709V"
## [1] "N70K"
## [1] "N70S"
## [1] "V714A"
## [1] "V714D"
## [1] "M715I"
## [1] "M715T"
## [1] "P717L"
## [1] "F718L"
## [1] "I723V"
## [1] "T724I"
## [1] "V728I"
## [1] "p.L729del"
## [1] "N730K"
## [1] "T731I"
## [1] "F733L"
## [1] "M734V"
## [1] "A735E"
## [1] "A735T"
## [1] "A735V"
## [1] "L736P"
## [1] "p.N740del"
## [1] "p.M741_T742delinsI "
## [1] "T742A"
## [1] "E746K"
## [1] "E747A"
## [1] "M748I"
## [1] "E74D"
## [1] "E74K"
## [1] "Q750R"
## [1] "V751F"
## [1] "V751I"
## [1] "G752R"
## [1] "G758E"
## [1] "I759V"
## [1] "E763D"
## [1] "E763K"
## [1] "M764K"
## [1] "M764R"
## [1] "L771V"
## [1] "D772N"
## [1] "P773S"
## [1] "Y774C"
## [1] "Y774D"
## [1] "p.Y776del"
## [1] "F777L"
## [1] "Q779K"
## [1] "G77R"
## [1] "N782T"
## [1] "I783T"
## [1] "F784L"
## [1] "D785N"
## [1] "I788V"
## [1] "V789A"
## [1] "V789I"
## [1] "L791F"
## [1] "L793F"
## [1] "G797V"
## [1] "P79A"
## [1] "P79H"
## [1] "P79S"
## [1] "P7L"
## [1] "P7R"
## [1] "R800C"
## [1] "R800H"
## [1] "R800L"
## [1] "M801V"
## [1] "p.801_803delMSN/insS"
## [1] "S805L"
## [1] "V806M"
## [1] "R808C"
## [1] "R808H"
## [1] "R808P"
## [1] "R811G"
## [1] "R811H"
## [1] "L812Q"
## [1] "R814Q"
## [1] "F816L"
## [1] "F816Y"
## [1] "K817E"
## [1] "W822C"
## [1] "P823T"
## [1] "N826D"
## [1] "L828V"
## [1] "D82E"
## [1] "G833R"
## [1] "N834D"
## [1] "S835A"
## [1] "S835L"
## [1] "V836M"
## [1] "L839P"
## [1] "N841K"
## [1] "T843A"
## [1] "L846R"
## [1] "I848F"
## [1] "D84G"
## [1] "D84N"
## [1] "V850M"
## [1] "F851L"
## [1] "V856L"
## [1] "G857D"
## [1] "M858L"
## [1] "P85S"
## [1] "P85T"
## [1] "S866L"
## [1] "S866P"
## [1] "E867K"
## [1] "E867Q"
## [1] "R869S"
## [1] "F86L"
## [1] "D872N"
## [1] "S873A"
## [1] "G874D"
## [1] "R878C"
## [1] "R878H"
## [1] "R878L"
## [1] "W879R"
## [1] "Y87C"
## [1] "H886P"
## [1] "H886Q"
## [1] "S88G"
## [1] "I890T"
## [1] "I891N"
## [1] "I891T"
## [1] "F892I"
## [1] "R893C"
## [1] "R893H"
## [1] "I894M"
## [1] "L895F"
## [1] "C896S"
## [1] "G897E"
## [1] "G897R"
## [1] "R8P"
## [1] "R8Q"
## [1] "R8W"
## [1] "E901K"
## [1] "S901L"
## [1] "S910L"
## [1] "G911E"
## [1] "Q912R"
## [1] "C915R"
## [1] "L917R"
## [1] "L917V"
## [1] "V922I"
## [1] "V924I"
## [1] "I925F"
## [1] "N927K"
## [1] "N927S"
## [1] "L928P"
## [1] "T92I"
## [1] "L935P"
## [1] "L939F"
## [1] "F93S"
## [1] "S940N"
## [1] "S941F"
## [1] "S941N"
## [1] "S943N"
## [1] "D945G"
## [1] "I94S"
## [1] "I94V"
## [1] "E952K"
## [1] "D953E"
## [1] "D953Y"
## [1] "R954G"
## [1] "M956I"
## [1] "M956T"
## [1] "L959P"
## [1] "V95I"
## [1] "V95L"
## [1] "Q960K"
## [1] "A964G"
## [1] "R965C"
## [1] "R965H"
## [1] "R965L"
## [1] "Q967R"
## [1] "G969C"
## [1] "G969S"
## [1] "R971C"
## [1] "R971H"
## [1] "K974D"
## [1] "R975Q"
## [1] "R975W"
## [1] "D979H"
## [1] "C981F"
## [1] "C982R"
## [1] "G983D"
## [1] "R986L"
## [1] "R986Q"
## [1] "R986W"
## [1] "R988Q"
## [1] "R988W"
## [1] "K991E"
## [1] "K991T"
## [1] "A993S"
## [1] "A993T"
## [1] "A994T"
## [1] "A994V"
## [1] "L995F"
## [1] "A996T"
## [1] "A997D"
## [1] "A997S"
## [1] "A997T"
## [1] "G999D"
## [1] "G9S"
## [1] "G9V"Step 6 Fit posterior mean \(\tilde\mu_i\) with clinical covariates
The clinical covariates include eaRate, blastpssm, provean_score, pph2_prob, ipeak, and penetrance density.
Here pattern mixture model is used.
Step 7 Use fitted penetrance and MSE from Step 6 to solve for \(\alpha\) and \(\beta\)
Adjust \(\alpha_{fi}\) and \(\beta_{fi}\) to \(\alpha_0\) and \(\beta_0\) if they were too small (<0.01)
Step 8 Calculate convergence criterion \(\delta\) to see if it’s small enough to end the loop
Calculate the convergence criterion \(\delta\) as follows, and if it’s greater than 0.1, then repeat Step 6 to Step 8. When \(\delta\) is smaller than 0.1, EM converges.
\[\delta=\sum_{i=1}^m |p^{em+}_i-p^{em-}_i| \]
where \(p_i^{em+}\) is the penetrance after one iteration of EM alogrithm for each variant, and \(p_i^{em-}\) is the penetrance from the previous iteration of EM alogrithm for each variant.
EM algorithm converges in 7 steps.
delta <- 10
d$eaRate <- as.numeric(d$eaRate)
d$blastpssm <- as.numeric(d$blastpssm)
d$ipeak <- suppressWarnings(as.numeric(d$ipeak))
covariates <- c("eaRate","blastpssm","provean_score","pph2_prob","ipeak","feat_dist_w")
# weight for weighted regression
d$weight_r <- 1/p_variance_w
d$alpha_1_w <- alpha_1_w
d$beta_1_w <- beta_1_w
d$p_mean_w <- p_mean_w
regression <- function(dv, pivs, nivs, data) {
# run a linear model with text arguments for dv and ivs
piv_string <- paste(pivs, collapse=" + ")
niv_string <- paste(nivs, collapse=" - ")
if(niv_string!="") iv_string <- paste(piv_string, " - ", niv_string, sep = "")
if(niv_string=="") iv_string <- paste(piv_string)
#print(iv_string)
regression_formula <- as.formula(paste(dv, iv_string, sep=" ~ "))
#print(regression_formula)
lm(regression_formula, data, weights = d$weight_r)
}
count <- 0
while(delta>0.1 & count < 10){
count <- count + 1
alpha_f <- NULL
beta_f <- NULL
for(i in 1:nrow(d)){
newdata = data.frame(var=d[i,"var"])
newdata[covariates] <- d[i,covariates]
model <- regression("p_mean_w", covariates,
colnames(newdata)[colSums(is.na(newdata))>0], d)
mean_f <- predict(model, newdata)
variance_f <- (predict(model, newdata,se.fit = T)$se.fit)^2
alpha <- solab(mean_f,variance_f)[1]
beta <- solab(mean_f,variance_f)[2]
if(alpha<0.01 | beta<0.01){
alpha_f[i]=alpha_0_w
beta_f[i]=beta_0_w
}else{
alpha_f[i]=alpha
beta_f[i]=beta
}
}
new_mean <- (alpha_f + d$brs1)/(alpha_f + beta_f + d$total_carriers)
delta <- sum(abs(new_mean-d$p_mean_w))
d$p_mean_w <- new_mean
}Scale the variance by a factor of 20
# update alpha and beta
# when tuning parameter is 20
mean <- alpha_f/(alpha_f+beta_f)
variance <- mean*(1-mean)
variance <- variance / 20
alpha_g <- (mean^2 * (1-mean) - variance * mean) /variance
beta_g <- alpha_g * (1 / mean - 1)
# save data
d$alpha_g <- alpha_g
d$beta_g <- beta_g
save(d,file="CheckPoint_20.RData")Results
Part 1 Coverage plots
Bootstrap and get the coverage rate
Use the observed penetrance from some variant as the TRUE penetrance for that variant, generate n binomial observations
Use the final EM algorithm posterior as the prior for Beta-Binomial, incorporate data from previous step, generate the posterior distribution, get 95% credible interval.
Check whether the interval cover the true penetrance from Step 1.
Repeat Step 1 to Step 3 N times to get the coverage rate.
Bootstrap function
BootsCoverage <- function(var,n=100,N=1000,true){
# var: variant name
# n: number of subjects in the new data
# N: number of Bootstrap
if(true=="empirical") true="p_empirical"
if(true=="observed") true="p_observed"
# extract the "true" penetrance
true.p <- d[d$var==var,true]
# generate binomial data
event <- rbinom(N,n,true.p)
# get the posterior credible interval
alpha <- d$alpha_g[which(d$var==var)]
beta <- d$beta_g[which(d$var==var)]
new.alpha <- alpha + event
new.beta <- beta + n - event
lb <- qbeta(0.025,new.alpha,new.beta)
ub <- qbeta(0.975,new.alpha,new.beta)
# change lb to floor of nearest 0.01
lb <- floor(lb*100)/100
return(sum(lb < true.p & ub > true.p)/N)
}Plot coverage
Empirical penetrance as the true penetrance
The coverage plot where the empirical penetrance is the true penetrance and one hundred new observations are added is shown below.
load("CheckPoint_20.RData")
########## when n = 100 ##########
results <- sapply(d$var,function(x) BootsCoverage(x,n=100,true="empirical") )
carriers.size <- ifelse(d$total_carriers <= 10,1,ifelse(d$total_carriers <= 100,2,ifelse(d$total_carriers <= 1000,3,ifelse(d$total_carriers <= 10^4,4,5) )))
new.data <- data.frame(Penetrance=d$p_empirical,Coverage=results,Number=log10(d$total_carriers))
ggplot(data=new.data,aes(x=Penetrance,y=Coverage))+geom_point(aes(size=Number,color=Number),shape=20)+geom_hline(yintercept = 0.95,color="red")+scale_x_continuous(limits=c(0,1))+scale_y_continuous(limits=c(0,1))+labs(x=" True penetrance under simulation", y="Coverage rate", size=TeX("$\\log_{10}$(Total number of carriers)"),color=TeX("$\\log_{10}$(Total number of carriers)"))+
theme(legend.position = "bottom",legend.box = 'vertical',legend.justification = 'left',legend.box.just = 'left',legend.title = element_text(size=8))+scale_colour_gradient(low = "dodgerblue", high = "black")
Observed penetrance as the true penetrance
The coverage plot where the observed penetrance is the true penetrance and one hundred new observations are added is shown below.
########## when n = 100 ##########
results <- sapply(d$var,function(x) BootsCoverage(x,n=100,true="observed") )
carriers.size <- ifelse(d$total_carriers <= 10,1,ifelse(d$total_carriers <= 100,2,ifelse(d$total_carriers <= 1000,3,ifelse(d$total_carriers <= 10^4,4,5) )))
new.data <- data.frame(Penetrance=d$p_observed,Coverage=results,Number=log10(d$total_carriers))
ggplot(data=new.data,aes(x=Penetrance,y=Coverage))+geom_point(aes(size=Number,color=Number))+geom_hline(yintercept = 0.95,color="red")+scale_x_continuous(limits=c(0,1))+scale_y_continuous(limits=c(0,1))+labs(x=" True penetrance under simulation", y="Coverage rate", size=TeX("$\\log_{10}$(Total number of carriers)"),color=TeX("$\\log_{10}$(Total number of carriers)"))+
theme(legend.position = "bottom",legend.box = 'vertical',legend.justification = 'left',legend.box.just = 'left',legend.title = element_text(size=8))+scale_colour_gradient(low = "dodgerblue", high = "black")
Part 2 Output Bland-Altman plots
EM prior and EM posterior
# EM prior and EM posterior
prior.mean <- d$alpha_g/(d$alpha_g+d$beta_g)
posterior.mean <- (d$alpha_g + d$brs1)/(d$alpha_g+d$beta_g+d$total_carriers)
ba.data <- data.frame(Difference=prior.mean-posterior.mean, Average=(prior.mean+posterior.mean)/2,Size=log10(d$total_carriers))
ggplot(data=ba.data,aes(x=Average,y=Difference))+geom_point(aes(size=Size,color=Size))+scale_x_continuous(limits=c(0,1))+scale_y_continuous(limits=c(-0.6,0.6),breaks=seq(-0.6,0.6,length=9))+labs(x="Average of EM prior penetrance and EM posterior penetrance",y="EM prior penetrance- EM posterior penetrance",size=TeX("$\\log_{10}$(Total number of carriers)"),color=TeX("$\\log_{10}$(Total number of carriers)"))+
theme(legend.position = "bottom",legend.box = 'vertical',legend.justification = 'left',legend.box.just = 'left',legend.title = element_text(size=10))+scale_colour_gradient(low = "dodgerblue", high = "black")
EM prior and empirical posterior
# EM prior and empirical posterior
prior.mean <- d$alpha_g/(d$alpha_g+d$beta_g)
posterior.mean <- (0.4492651 + d$brs1)/(0.4492651+2.727777+d$total_carriers)
ba.data <- data.frame(Difference=prior.mean-posterior.mean, Average=(prior.mean+posterior.mean)/2,Size=log10(d$total_carriers))
ggplot(data=ba.data,aes(x=Average,y=Difference))+geom_point(aes(size=Size,color=Size))+scale_x_continuous(limits=c(0,1))+scale_y_continuous(limits=c(-0.6,0.6),breaks=seq(-0.6,0.6,length=9))+labs(x="Average of EM prior penetrance and empirical posterior penetrance",y="EM prior penetrance- Empirical posterior penetrance",size=TeX("$\\log_{10}$(Total number of carriers)"),color=TeX("$\\log_{10}$(Total number of carriers)"))+
theme(legend.position = "bottom",legend.box = 'vertical',legend.justification = 'left',legend.box.just = 'left',legend.title = element_text(size=10))+scale_colour_gradient(low = "dodgerblue", high = "black")+guides(color = guide_colorbar(order = 1),
size = guide_legend(order = 2))
EM posterior and empirical posterior
# EM posterior and empirical posterior
prior.mean <- (d$alpha_g + d$brs1)/(d$alpha_g+d$beta_g+d$total_carriers)
posterior.mean <- (0.4492651 + d$brs1)/(0.4492651+2.727777+d$total_carriers)
ba.data <- data.frame(Difference=prior.mean-posterior.mean, Average=(prior.mean+posterior.mean)/2,Size=log10(d$total_carriers))
ggplot(data=ba.data,aes(x=Average,y=Difference))+geom_point(aes(size=Size,color=Size))+scale_x_continuous(limits=c(0,1))+scale_y_continuous(limits=c(-0.6,0.6),breaks=seq(-0.6,0.6,length=9))+labs(x="Average of EM posterior penetrance and empirical posterior penetrance",y="EM posterior penetrance- Empirical posterior penetrance",size=TeX("$\\log_{10}$(Total number of carriers)"),color=TeX("$\\log_{10}$(Total number of carriers)"))+
theme(legend.position = "bottom",legend.box = 'vertical',legend.justification = 'left',legend.box.just = 'left',legend.title = element_text(size=10))+scale_colour_gradient(low = "dodgerblue", high = "black")
Part 3 Forest plots
The code to produce the forest plots is shown in the code chunk.
rm(d)
load("CheckPoint_20.RData")
new_mean <- (d$alpha_g + d$brs1)/(d$alpha_g+d$beta_g+d$total_carriers)
alt.mean <- (0.4492651 + d$brs1)/(0.4492651+2.727777+d$total_carriers)
lower <- qbeta(0.025,d$alpha_g,d$beta_g)
higher <- qbeta(0.975,d$alpha_g,d$beta_g)
forest.data <- data.frame(variant = d$var, em.posterior = new_mean, empirical.posterior=alt.mean,
lower=lower, higher=higher,total.carriers=d$total_carriers)
forest.data <- forest.data[order(forest.data$variant,forest.data$total.carriers),]
head(forest.data)
plotf <- function(a,b){
png( paste("Plots/Forest Plots/", a, "-",b,"pics.png",sep=""),res=300,height=10,width=10,units="in")
forestplot(paste("Variant: ",forest.data[a:b,"variant"],", total carriers:",forest.data$total.carriers[a:b]),
fn.ci_norm = fpDrawNormalCI,
boxsize = .25, # We set the box size to better visualize the type
line.margin = .1, # We need to add this to avoid crowding
mean = cbind(forest.data[a:b,"em.posterior"]),
lower = cbind(forest.data[a:b,"lower"]),
upper = cbind(forest.data[a:b,"higher"]),
xticks=seq(0,1,0.1),
col=fpColors(box="blue",line="red"),
xlab="Number line")
dev.off()
}
sapply(0:27*50+1,function(x) plotf(x,x+49) )
plotf(1401,1439)The forests plots are pasted below, where blue boxes are EM posteriors, and red segements are EM priors.
Part 4 Correlation Coefficients and Plots
The code to produce the weighted correlations is shown in the code chunk. The size of the circles in the presented figures are scaled by the Log10(Total Number of Carriers).
## wCorr v1.9.1load("CheckPoint_20.RData")
d<-d[,!(names(d) %in% "aaneigh_rank_inc_loops")]
plt.disease <- function(d, func, funcName, disName, xl=c(0,200), yl=c(0,100), variant="NA", sp=0.7){
par(cex=1, bty='l', lwd=2)
if (disName=="BrS_penetranceBayesian"){dis="BrS1"; pcolor="gray"}
else if(disName=="BrS_penetranceBayesian_initial"){dis="BrS1"; pcolor="gray"}
else if (disName=="all_penetranceBayesian"){dis="BrS1 and LQT3"; pcolor="gray"}
plot(d[,func], (d[,disName]), pch=21,
bg=pcolor, cex=log(d[,"total_carriers"])+1, lwd=1.5, axes=FALSE,
cex.lab=1.5, ylab=paste("Penetrance (%",dis,")",sep=""), xlab = funcName,
ylim = yl, xlim = xl)
axis(side=1,lwd=3, cex.axis=1.5)
axis(side=2,lwd=3, cex.axis=1.5)
abline(a=0,b=1)
if (variant!="NA"){
points(d[d$var==variant,func], 100*(d[d$var==variant,disName]), pch=21, cex=log(d[d$var==variant,"total_carriers"])+1, lwd=1.5, col="black", bg="red")
}
}
d$BrS_penetranceBayesian<-(d$alpha_g+d$brs1)/(d$alpha_g+d$beta_g+d$total_carriers)
filt<-0
fglm<-d[d$total_carriers>filt & !is.na(d$blastpssm) & !is.na(d$ipeak) & !is.na(d$feat_dist_w), ]
model <- lm(BrS_penetranceBayesian ~ ipeak, fglm, weights = 1/fglm$p_variance_w)
fglm$pred.ipeak<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.ipeak", "Predicted Penetrance", "BrS_penetranceBayesian", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Peak Current Covariate Alone"## [1] 0.3519097mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.3519097## 2.5% 97.5%
## 0.2436698 0.4620624## [1] 2.00000 17.79004model <- lm(BrS_penetranceBayesian ~ feat_dist_w, fglm, weights = 1/fglm$p_variance_w)
fglm$pred.feat_dist<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.feat_dist", "Predicted Penetrance", "BrS_penetranceBayesian", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Penetrance Density Covariate Alone"## [1] 0.6617035mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.6617035## 2.5% 97.5%
## 0.5228591 0.7578649## [1] 2.0000 -128.4841model <- lm(BrS_penetranceBayesian~ipeak+feat_dist_w,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.ipeak.dist<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.ipeak.dist", "Predicted Penetrance", "BrS_penetranceBayesian", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Peak Current and Penetrance Density Covariates"## [1] 0.7562233mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.7562233## 2.5% 97.5%
## 0.6540318 0.8305172## [1] 3.0000 -200.2098model <- lm(BrS_penetranceBayesian~ipeak+feat_dist_w+eaRate+blastpssm+provean_score+pph2_prob,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.all.comps.ipeak<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.all.comps.ipeak", "Predicted Penetrance", "BrS_penetranceBayesian", xl=c(0,0.9),yl=c(0,0.9))
## [1] "In Silico, Peak Current, and Penetrance Density Covariates"## [1] 0.7806726mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.7806726## 2.5% 97.5%
## 0.6937503 0.8544453## [1] 7.0000 -215.9894model <- lm(BrS_penetranceBayesian~eaRate+blastpssm+provean_score+pph2_prob,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.all.comps<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.all.comps", "Predicted Penetrance", "BrS_penetranceBayesian", xl=c(0,0.9),yl=c(0,0.9))
## [1] "In Silico Covariates Alone"## [1] 0.1942237mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.1942237## 2.5% 97.5%
## 0.1244528 0.2776599## [1] 5.00000 72.78969# Empirical
d$BrS_penetranceBayesian_initial<-(alpha_0_w+d$brs1)/(alpha_0_w+beta_0_w+d$total_carriers)
filt<-0
fglm<-d[d$total_carriers>filt & !is.na(d$blastpssm) & !is.na(d$ipeak) & !is.na(d$feat_dist_w), ]
model <- lm(BrS_penetranceBayesian_initial ~ ipeak, fglm, weights = 1/fglm$p_variance_w)
fglm$pred.ipeak<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.ipeak", "Predicted Penetrance", "BrS_penetranceBayesian_initial", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Peak Current Covariate Alone"## [1] 0.2226512mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian_initial,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian_initial[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.2226512## 2.5% 97.5%
## 0.1296066 0.3266152## [1] 2.0000 155.3884model <- lm(BrS_penetranceBayesian_initial ~ feat_dist_w, fglm, weights = 1/fglm$p_variance_w)
fglm$pred.feat_dist<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.feat_dist", "Predicted Penetrance", "BrS_penetranceBayesian_initial", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Penetrance Density Covariate Alone"## [1] 0.3513447mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian_initial,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian_initial[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.3513447## 2.5% 97.5%
## 0.1826644 0.4897752## [1] 2.0000 114.6661model <- lm(BrS_penetranceBayesian_initial~ipeak+feat_dist_w,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.ipeak.dist<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.ipeak.dist", "Predicted Penetrance", "BrS_penetranceBayesian_initial", xl=c(0,0.9),yl=c(0,0.9))
## [1] "Peak Current and Penetrance Density Covariates"## [1] 0.4225387mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian_initial,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian_initial[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.4225387## 2.5% 97.5%
## 0.2794866 0.5603973## [1] 3.00000 90.50763model <- lm(BrS_penetranceBayesian_initial~ipeak+feat_dist_w+eaRate+blastpssm+provean_score+pph2_prob,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.all.comps.ipeak<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.all.comps.ipeak", "Predicted Penetrance", "BrS_penetranceBayesian_initial", xl=c(0,0.9),yl=c(0,0.9))
## [1] "In Silico, Peak Current, and Penetrance Density Covariates"## [1] 0.436451mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian_initial,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian_initial[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.436451## 2.5% 97.5%
## 0.2783023 0.5644973## [1] 7.00000 93.02051model <- lm(BrS_penetranceBayesian_initial~eaRate+blastpssm+provean_score+pph2_prob,fglm, weights = 1/fglm$p_variance_w)
fglm$pred.all.comps<-predict(model, fglm, type="response")
plt.disease(fglm, "pred.all.comps", "Predicted Penetrance", "BrS_penetranceBayesian_initial", xl=c(0,0.9),yl=c(0,0.9))
## [1] "In Silico Covariates Alone"## [1] 0.1174166mod<-data.frame(model$fitted.values,model$model$BrS_penetranceBayesian_initial,model$weights)
foo <- boot(mod, function(data,indices)
weightedCorr(mod$model.fitted.values[indices],mod$model.model.BrS_penetranceBayesian_initial[indices], method="pearson", weights = mod$model.weights[indices])^2, R=1000)
foo$t0## [1] 0.1174166## 2.5% 97.5%
## 0.05564675 0.19489625## [1] 5.0000 189.9553